Random restoration projects

Summary plots

# total random sites to create
tot <- nrow(restdat)

id <- stri_rand_strings(tot, 5)
dts <- range(restdat$date)

ext <- bbox(tbpoly)

lon <- rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
lat <- rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
tmp <- SpatialPoints(cbind(lon, lat), 
                     proj4string = crs(tbpoly)
) %>% 
  .[tbpoly, ] %>% 
  .[sample(1:nrow(.@coords), tot, replace = F), ]

restdat_rnd <- tibble(id) %>% 
  mutate(
    date = rnorm(tot, mean(restdat$date), sd(restdat$date)),
    date = round(date, 0),
    top = sample(c('hab', 'wtr'), tot, replace = T)    
  )

reststat_rnd <- tibble(id) %>% 
  mutate(
    lat = tmp$lat, 
    lon = tmp$lon
  )

resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

# map by restoration type
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)

# map by date
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)

# barplot of date counts
toplo <- restall_rnd %>% 
  group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) + 
  geom_bar() + 
  coord_flip() + 
  theme_bw() + 
  theme(
    axis.title.y = element_blank()
  ) +
  scale_y_discrete(expand = c(0, 0))

Distance to restoration sites

wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')

head(wqmtch_rnd)

## # A tibble: 6 x 5
##    stat resgrp   rnk    id     dist
##   <int>  <chr> <int> <chr>    <dbl>
## 1    47    hab     1 IYJz0 2960.310
## 2    47    hab     2 3OOJT 3394.566
## 3    47    hab     3 Jp5pr 3549.903
## 4    47    hab     4 GD6tr 4737.298
## 5    47    hab     5 nkm7v 4759.014
## 6    47    hab     6 Vw1c6 5440.625

Closest

## 
# plots

# combine lat/lon for the plot
toplo <- wqmtch_rnd %>% 
  left_join(wqstat, by = 'stat') %>% 
  left_join(reststat_rnd, by = 'id') %>% 
  rename(
    `Restoration\ngroup` = resgrp,
    `Distance (dd)` = dist
  )
    
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")

# base map
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
  geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)

# closest
toplo1 <- filter(toplo, rnk %in% 1)

pbase + 
  geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest twenty percent

# closest five percent
fvper <- max(toplo$rnk) %>% 
  `*`(0.2) %>% 
  ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))

pbase + 
  geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest all combinations

# closest all combo
toplo3 <- toplo

pbase + 
  geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Summarizing effects of restoration projects

Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.

salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)

## # A tibble: 6 x 3
## # Groups:   stat [2]
##    stat     cmb     cval
##   <int>   <chr>    <dbl>
## 1     6 hab_aft 25.09977
## 2     6 hab_bef 24.25426
## 3     6 wtr_aft 24.05831
## 4     6 wtr_bef 24.86349
## 5     7 hab_aft 24.82157
## 6     7 hab_bef 24.98083

head(salchg_rnd)

## # A tibble: 6 x 4
##    stat     hab     wtr     cval
##   <int>  <fctr>  <fctr>    <dbl>
## 1     6 hab_aft wtr_aft 24.57904
## 2     6 hab_aft wtr_bef 24.98163
## 3     6 hab_bef wtr_aft 24.15629
## 4     6 hab_bef wtr_bef 24.55887
## 5     7 hab_aft wtr_aft 24.62957
## 6     7 hab_aft wtr_bef 25.25801

Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.

wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)

## # A tibble: 4 x 5
##       hab     wtr              data       crv                    prd
##    <fctr>  <fctr>            <list>    <list>                 <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>

Discretization of salinity conditional probability distributions:

salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd

## # A tibble: 8 x 5
##       hab     wtr      qts       brk  clev
##    <fctr>  <fctr>    <dbl>     <dbl> <dbl>
## 1 hab_aft wtr_aft 26.52550 0.4091514     1
## 2 hab_aft wtr_aft 30.03970 0.8643254     2
## 3 hab_aft wtr_bef 25.88986 0.3671466     1
## 4 hab_aft wtr_bef 29.78396 0.8493238     2
## 5 hab_bef wtr_aft 26.85391 0.4598276     1
## 6 hab_bef wtr_aft 30.24440 0.8849303     2
## 7 hab_bef wtr_bef 26.21855 0.4170020     1
## 8 hab_bef wtr_bef 29.98880 0.8717005     2

A plot showing the breaks:

toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>% 
  unnest
ggplot(toplo, aes(x = cval, y = cumest)) + 
  geom_line() + 
  geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr) +
  theme_bw()

Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.

# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')

# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>% 
  group_by(hab, wtr) %>% 
  nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>% 
  rename(
    salev = cval.y, 
    cval = cval.x
  ) %>% 
  group_by(hab, wtr) %>% 
  nest %>% 
  left_join(salbrk_rnd, by = c('hab', 'wtr')) %>% 
  mutate(
    sallev = pmap(list(data, levs), function(data, levs){
      # browser()
      out <- data %>% 
        mutate(
          saval = salev,
          salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          salev = as.character(salev)
        )
      
      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')

chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')

chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>% 
  print(n = nrow(.))

## # A tibble: 24 x 6
##        hab     wtr salev       qts       brk  clev
##     <fctr>  <fctr> <chr>     <dbl>     <dbl> <dbl>
##  1 hab_aft wtr_aft    lo 10.914175 0.6595735     1
##  2 hab_aft wtr_aft    lo 15.669895 0.9759047     2
##  3 hab_aft wtr_aft    md  6.878474 0.5760369     1
##  4 hab_aft wtr_aft    md  9.498920 0.9468640     2
##  5 hab_aft wtr_aft    hi  3.584727 0.3016759     1
##  6 hab_aft wtr_aft    hi  4.363590 0.7610040     2
##  7 hab_aft wtr_bef    lo  9.092627 0.4355175     1
##  8 hab_aft wtr_bef    lo 12.363819 0.8815996     2
##  9 hab_aft wtr_bef    md  6.103452 0.3793212     1
## 10 hab_aft wtr_bef    md  8.129952 0.8387278     2
## 11 hab_aft wtr_bef    hi  3.473219 0.2582490     1
## 12 hab_aft wtr_bef    hi  4.049840 0.7176782     2
## 13 hab_bef wtr_aft    lo 11.089097 0.6925214     1
## 14 hab_bef wtr_aft    lo 16.592253 0.9863017     2
## 15 hab_bef wtr_aft    md  5.667113 0.3303518     1
## 16 hab_bef wtr_aft    md  6.648293 0.7580916     2
## 17 hab_bef wtr_aft    hi  4.077344 0.4627592     1
## 18 hab_bef wtr_aft    hi  5.229716 0.8806475     2
## 19 hab_bef wtr_bef    lo  9.858977 0.5375839     1
## 20 hab_bef wtr_bef    lo 13.586304 0.9302304     2
## 21 hab_bef wtr_bef    md  5.938513 0.4033686     1
## 22 hab_bef wtr_bef    md  7.217921 0.8247996     2
## 23 hab_bef wtr_bef    hi  3.616690 0.3122743     1
## 24 hab_bef wtr_bef    hi  4.315727 0.7382746     2

Final combinations long format:

chlbar_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  mutate(
    data = map(data, function(x){
      
      brk <- x$brk
      out <- data.frame(
        lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
      )
      
      return(out)
      
    })
  ) %>% 
  unnest %>% 
  gather('chllev', 'chlval', lo:hi) %>% 
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi')),
    chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
  )
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')

chlbar_rnd %>% 
  print(n = nrow(.))

## # A tibble: 36 x 5
##        hab     wtr  salev chllev     chlval
##     <fctr>  <fctr> <fctr> <fctr>      <dbl>
##  1 hab_aft wtr_aft     lo     lo 0.65957353
##  2 hab_aft wtr_aft     md     lo 0.57603686
##  3 hab_aft wtr_aft     hi     lo 0.30167585
##  4 hab_aft wtr_bef     lo     lo 0.43551747
##  5 hab_aft wtr_bef     md     lo 0.37932121
##  6 hab_aft wtr_bef     hi     lo 0.25824904
##  7 hab_bef wtr_aft     lo     lo 0.69252135
##  8 hab_bef wtr_aft     md     lo 0.33035181
##  9 hab_bef wtr_aft     hi     lo 0.46275919
## 10 hab_bef wtr_bef     lo     lo 0.53758390
## 11 hab_bef wtr_bef     md     lo 0.40336858
## 12 hab_bef wtr_bef     hi     lo 0.31227425
## 13 hab_aft wtr_aft     lo     md 0.31633115
## 14 hab_aft wtr_aft     md     md 0.37082709
## 15 hab_aft wtr_aft     hi     md 0.45932817
## 16 hab_aft wtr_bef     lo     md 0.44608217
## 17 hab_aft wtr_bef     md     md 0.45940662
## 18 hab_aft wtr_bef     hi     md 0.45942914
## 19 hab_bef wtr_aft     lo     md 0.29378034
## 20 hab_bef wtr_aft     md     md 0.42773980
## 21 hab_bef wtr_aft     hi     md 0.41788835
## 22 hab_bef wtr_bef     lo     md 0.39264652
## 23 hab_bef wtr_bef     md     md 0.42143103
## 24 hab_bef wtr_bef     hi     md 0.42600033
## 25 hab_aft wtr_aft     lo     hi 0.02409532
## 26 hab_aft wtr_aft     md     hi 0.05313605
## 27 hab_aft wtr_aft     hi     hi 0.23899597
## 28 hab_aft wtr_bef     lo     hi 0.11840036
## 29 hab_aft wtr_bef     md     hi 0.16127217
## 30 hab_aft wtr_bef     hi     hi 0.28232182
## 31 hab_bef wtr_aft     lo     hi 0.01369831
## 32 hab_bef wtr_aft     md     hi 0.24190839
## 33 hab_bef wtr_aft     hi     hi 0.11935246
## 34 hab_bef wtr_bef     lo     hi 0.06976958
## 35 hab_bef wtr_bef     md     hi 0.17520039
## 36 hab_bef wtr_bef     hi     hi 0.26172541

Discretesize chlorophyll data, all stations:

# discretize all chl data by breaks 
chlbrk_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>% 
  mutate(
    lev = pmap(list(data, levs), function(data, levs){
      
      out <- data %>% 
        mutate(
          lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          lev = as.character(lev)
        )

      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest %>% 
  rename(
    chlev = lev, 
    chval = cval
    )

save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')

A bar plot of splits:

ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  facet_grid(hab ~ wtr) +
  theme_bw()

A plot showing the breaks:

toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>% 
  unnest %>%
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi'))
  )
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) + 
  geom_line() + 
  geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr, scales = 'free_x') +
  theme_bw()